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In the coarse grained Brownian Dynamics simulation method the many solvent molecules are replaced by 
random thermal kicks and an effective friction acting on the particles of interest. For Brownian Dynamics the 
friction has to be so strong that the particles' velocities are damped much faster than the duration of an inte- 
gration timestep. Here we show that this conceptual limit can be dropped with an analytic integration of the 
equations of damped motion. In the resulting Langevin integration scheme our recently proposed approximate 
form of the hydrodynamic interactions between the particles can be incorparated conveniently, leading to a fast 
multi-particle propagation scheme, which captures more of the short-time and short-range solvent effects than 
standard BD. Comparing the dynamics of a bead-spring model of a short peptide, we recommend to run simu- 
lations of small biological molecules with the Langevin type finite damping and to include the hydrodynamic 
interactions. 
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I. INTRODUCTION 



For the simulation of microscopic particles in a solvent, 
Brownian Dynamics (BD) has become a work horse tech- 
nique, which is routinelvianplied to study, e.g., the pronartigs 
of colloidal suspensionaJJU^ the dynamicS|«f,p«lymersLLl]U, 
or the association of biological proteinsLUJUU. Einstein's 
seminal explanation of thejjandom thermal motion of micro- 
scopic objects in a solventLJ, which had been discovered ear- 
lier by Robert Brown, is based on the assumption that the sol- 
vent molecules are so much smaller than the interesting ob- 
jects that they can be replaced by a continuous solvent, which 
then modifies the motion of the much larger Brownian par- 
ticles in three major ways. First, the kinetic energy of the 
large particles is dissipated via the viscosity of the solvent. 
Second, the thermal motion of the solvent molecules is re- 
placed by random kicks to the observed particles. These two 
contributions effectively form a thermostat. Third, the direct 
interactions like Coulomb or short range van-der-Waals inter- 
actions are modified by the solvent. With these three modifi- 
cations, the problem of solving Newton's equations of motion 
for very many particles with well-defined position-dependent 
interaction potentials is replaced by solving a few-body prob- 
lem with more complex interactions, a velocity-dependent 
dissipation, and an additional noise term. This level of ap- 
proximation is often referred to as Langevin Dynamics (LD), 
whereas Brownian Dynamics involves the additionaLapproxi- 
mation that only "long" time intervals are consideredU. Then, 
due to the damping of the particle velocities, all information 
about the momentum is lost between two subsequent observa- 
tions. This is called the overdamped regime. 

Obviously, the central criterion whether BD is an adequate 
simulation method, is that the solvent molecules are indeed 
much faster — or that the Brownian particles are much larger 
and thus slower than a water molecule with its effective diam- 
eter of about 1 .5 A. For colloidal particles with diameters on 
the micron scale this separation of scales is a good approxi- 
mation, but for the much smaller biological proteins with di- 
ameters of a few nanometers it is already questionable. This 



becomes even more pronounced, when this coarse grained 
method is applied to the internal dynamics of proteins. 

In addition to the thermostat effect, the solvent also medi- 
ates a mechanical coupling between the observable particles, 
the so called hydrodynamic interactions (HI). The resulting 
correlations in the particle velocities may have striking^ffects 
on the folding behavior of proteins as shown recentlyU. 

Here, we investigate the effects of the approximations out- 
lined above. For a bead-spring representation of a short pep- 
tide we compare plain BD simulations to BD simulations with 
HI included and also to LD simulations both with and without 
HI to see how the simulated dynamics of very small parti- 
cles is influenced by the non-negligible relaxation times and 
how the hydrodynamic coupling between the beads affects the 
overall dynamics. In the next section we introduce our imple- 
mentation of the LD and BD algorithms and how hydrody- 
namic interactions can be incorporated efficiently. There we 
also explain the coarse graining procedure used to set up the 
bead-spring model of the peptide. The results then show how 
the different approximations influence the dynamic behavior 
of the model peptide. 



II. METHODS 
A. Langevin and Brownian Propagation 

For the following it is convenient to formulate the Brownian 
and the Langevin propagation algorithms with effective forces 
instead of directly using displacempats as in the original BD 
scheme of Ermak and McCammonU. For this, we start from 
Newton's equations of motion for the complete system includ- 
ing all solvent molecules. Obviously, the particle masses rrii 
can be taken as constant such that the change of the momenta 
Pi can be expressed as a change of the velocities Vi due to the 
sums Fi of all pairwise position-dependent forces. 
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The dimensionality of this system of equations can be greatly 
reduced by assuming that the solvent molecules are very 
small. Then their masses are very small and they move so fast 
that their effect on the slower degrees of freedom of the inter- 
esting particles can be treated as a mean-field heat bath coa- 
sisting of a dissipative friction term plus a random driving|13. 
For a spherical particle, the phenomenological friction con- 
stant 7 can be calculated via the Stokes-Einstein relation from 
the hydrodynamic radius a of the particle aruLthe macroscopic 
bulk viscosity rj of the solvent as 7 = BTTTyaLj. Due to the ran- 
dom nature of the driving kicks from the solvent molecules, 
these kicks are characterized via the statistical measures of a 
vanishing average and a finite covariance of the resulting dis- 
placements Ri{dt): 

{R,) - and {R,Rk) = 2D,kSt (2) 

For the diagonal terms with i ^ k the diffusion coefficient 
Dii = ksT/^i gives the ratio between the solvent induced 
thermal energy fc^T and the also solvent mediated dissipa- 
tion. For the off-diagonal terms, Dik describes the hydrody- 
namic coupling due to the displaced solvent. 

To determine the random displacements Ri from the co- 
variance (|2|), one^ssentially has to take the square root of the 
diffusion matrixLl Due to the high numerical cost associated 
with this step, the hydrodynamic coupling was often neglected 
in BD simulations by setting the off-diagonal entries Di^. ~ 0. 
It was only recently that we could show how this expensive 
matrix factorization can be approximated by a sum of two- 
body contributions and the correlated random-displacements 
can be expressed via effective random forcesU f^^^, leading 
to the same C'(iV^) runtime scaling as for the pairwise direct 
interactions Fn- . 

To arrive at a formulation of the equation of motion with 
the implicit solvent which is convenient for numerical propa- 
gation, we start from ([ij) with the total force F = F^xt + f, 
which is the sum of the external forces F^xt and the random 
kicks /: 

'i = -iF~,v) (3) 
at m 

For simplicity, we omit the coordinate index i for the fol- 
lowing. This equation of motion can be integrated analyt- 
ically over one timestep of length At, which must be so 
small that the_configuration dependent forces remain essen- 
tially constantU. Due to the linear Stokes friction, we can use 
the average /(Ai) of the random kicks. With the initial ve- 
locity vq, v{At) at the end of the timestep is then 

F f F\ 7At 

v{At) = - + Uo - - e-— (4) 
7 V 7/ 

This equation can be integrated once more to give the dis- 
placement Ax{At) during At: 

A.(AO = -At-I^f^- J(l-e-^) (5) 
7 7 V7 / ^ ^ 

From these equations the conventional BD propagation al- 
gorithm is obtained by going to the overdamped regime where 




FIG. 1: Range of admissible timesteps At between the inaccurate 
and the unphysical regimes of a BD simulation vs. particle radius 
a. The lower green line denotes the conceptual limit according to 
At ^ Trei oc . The upper red limit indicates that the maximal 
timestep for a numerically stable propagation increases only oc a. 
Consequently, for any At there is a minimal particle size amin, for 
which the propagation becomes unstable, and a maximal Umax, for 
which the assumption of overdamped motion is not valid anymore. 
For further explanations see text. 

At ^ Trei, i.e., only timesteps At are considered that are 
much longer than the velocity relaxation time Trei ~ m/7. 
Then vq can be neglected and the displacement simplifies to 

Ax(At) = -At= -^FAt, (6) 
7 ksT 

which can also be expressed with the diffusion constant Dq = 
kgT/j of the particle. 

For the rotational motion, analogous equations are used 
where position x and velocity v are replaced by a rotation an- 
gle and an angular velocity, and torques instead of the forces 
F act on moments of inertia. 

The Langevin equations and look much more com- 
plex than the Brownian propagator (j^), but for practical appli- 
cations the main effort is to determine the forces and torques 
acting on each of the particles. Thus, keeping track of the ve- 
locity only incurs a negligibly small additional overhead for 
the benefit that one has not to make sure that the conceptual 
constraint of BD, At ^ Tr^i, is fulfilled. It is actually quite 
easy to find simulation scenarios where these BD require- 
ments cannot be fulfilled for all particles simultaneously. This 
can be seen from the following estimate of the two limits for 
At. The first is that At ^ Trei ~ m/7. For typical colloidal 
particles or biological polymers, the mass m = ^ pa'^ is pro- 
portional to the volume and thus scales oc with the radius 
a. From the Stokes-Einstein relation we get 7 = Girrja oc a, 
thus Trei — | incrcascs quadratically with a. On the other 
hand. At has to be small enough for numerical reasons to keep 
the typical displacements Aa;(At) <^ Axmax, where Axmax 
is a typical potential well extensions or the dimensions of the 
smallest particles. According to (||) this upper limit scales 
oc a. In a simulation where all particles have the same size, 
usually a At can be found, which fulfills both requirements. 
However, when particles of different sizes are considered, a 
timestep which allows for a stable propagation of the small 
particles may be much shorter than T^ei for the largest par- 
ticles, while a At that ensures the overdamped regime for 
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the larger particles most likely leads to numerical instabilities. 
This problem is sketched in figure |lj where the green line de- 
notes the lower limit from the overdamped regime and the red 
line the upper numerical stability margin. The dotted area be- 
tween these two limits is the area of permissible timesteps. 
For a fixed a there is a certain region [At„im, lS.t,nax\ for At 
to choose from, but when particles of different sizes are used 
in the same simulation, then for a given timestep At there is 
only a range of particle sizes [omm, cLmax] which can be prop- 
agated reliably within the assumption of Brownian Dynamics. 

A few numbers highlight this problem. With a typical pro- 
tein density of p = 1.1/9^2 and p/rj = 1.23 ps/nm^ we find 
that a protein of a = 5 nm radius has Trei = 6.8 ps. The small 
electron carrier cytochrome C2 with a = 1.67 nm has only 
Trei = 0.76 ps and an ion with an effective hydrodynamic ra- 
dius of 0.2 nm looses its velocity already within Trei ~ 0.01 
ps. In a simulation with cytochrome C2 and ions, based on 
practical experience, a timestep of At = 1 ps still yields a sta- 
ble propagation of the fast ions. While the ions are clearly in 
the overdamped regime, this assumption is at least question- 
able for the cytochrome C2 . When one wants to simulate the 
association of cytochrome C2 to larger proteins in the presence 
of explicitly modelled ions. At = 1 ps is clearly too short for 
a theoretically vahd overdamped BD description. 

We note that for the analytical integration above we re- 
quired that the interparticle forces be constant during one 
timestep. When the underlying interaction potentials are not 
approximated by constant slopes but one order further by har- 
monic potentials, the analytical solutions for the motion of 
a damped harmonic oscillator can be used for the propaga- 
tor. Though this is not really practical, it allows to distin- 
guish between the strongly damped creeping and the weakly 
damped ballistic regimes based on the ratio of potential curva- 
ture and relaxation time. With this criterion and the parameter 
values above, only for very short ranged potentials like van- 
der-Waals interactions the dynamics may come close to the 
weakly damped balUstic and thus oscillatory regime. 



B. Including Hydrodynamic Interactions 

When one wants to include HI into a Langevin propaga- 
tion as the one from equations andp(E|), then there is a 
conceptual|-(Mficulty. The usual OseenLjor Rotne-Prager- 
Yamakawal— U (RPY) hydrodynamic tensors are built with the 
help of the Faxen theoremU from the stationary flow fields 
around spheres moving with constant velocities. In an over- 
damped BD simulation one can at least assume that the ve- 
locities are constant during one timestep and then treat each 
timestep as a configuration with (temporarily) stationary ve- 
locities mimicking these creeping flow conditions. Due to the 
linear Stokes friction, the velocity does not even occur in the 
BD propagator (H) and the (constant) forces can directly be 
converted into displacements. 

With hydrodynamic intepajCtions the BD propagator ^ for 
the ith coordinate becomesllJ 

Ax,{At) = V ^^^At + R,iAt), (7) 



fc 



which can bejrewritten with a hydrodynamically corrected ef- 
fective forceO Fl;^'^ = J2k %^Fk- 
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(8) 



With the RPY HI tensor, the term taking care of local vari- 
ations of the diffusion coefficients vanishes and is omit- 
ted here (see, e.g., the original derivation r-by Ermak and 
McCammonLj). Recently, we could showU how the ex- 
pensive calculation of the correlated random displacements 
Ri{At) can be approximated efficiently with an ansatz that 
converts the uncorrelated random forces fi into hydrodynam- 
ically corrected effective random forces 
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With f3ii — 1, the normalization factors can be determined 
from 



and the quadratic equation 
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1 - y/l ~ [(TV - 1)62 _ (jV _ 2)g] 

(iV - l)e2 - (TV - 2)e 



(11) 



where e — {Dik/ Da) is the average normalized off-diagonal 
entry of the diffusion matrix. Then, the displacements Axi 
can be calculated efficiently with an overall 0{N'^) runtime 
scaling from the hydrodynamically corrected external and ran- 
dom forces and the diagonal entries of the diffusion matrix as 



Axi{At) 



Du At . 
ksT 



(Ff' + ft^n 



(12) 



We now show how this idea can also be used in an LD 
description, where the velocities are not constant during a 
timestep but change from the initial uo to the final w(At), 
which is then the new wq for the next step. For this, we need 
an effective force F which is constant during one timestep 
and leads to the same total displacement as with the correct 
propagator With A.T(At) = FAt/^ we find from equation 
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The new vq for the next timestep is then 
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(13) 
(14) 



In the actual implementation equations (jlj) and (jlj) occur 
twice, once for the external forces Fi and once for the random 
forces fi. From these damping-corrected uncorrelated forces 
the hydrodynamically correlated forces F[-^-^ and f^-^^ are 
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FIG. 2: X-component of the trajectory of the peptide CM under 
the same sequence of random kicks with the different propagation 
schemes. The inclusion of HI accelerates the CM diffusion and the 
finite damping of LD smoothes and slightly delays the peptide's re- 
sponse to the fast thermal kicks without changing the long-time tra- 
jectory. 



then computed as explained above and used to propagate the 
particles according to equation ([I^).^ith this recipe, our re- 
cent efficient truncated expansion HlO can be combined with 
a Langevin type propagation of Brownian particles, which is 
conceptually valid even for arbitrarily small timesteps with 
only negligibly more effort. Most of the prefactors in equa- 
tions (^3]) and ( p^ need only be evaluated once during the 
simulation setup. 



III. RESULTS 

As a testcase for the propagation schemes outlined above 
we used a small eleven-residue peptide, which had been used 
in our group in a recent molecular dynamics (MD) study of the 
effects of interfacial water-layers during peptide-protein asso- 
ciation at an SH3 domainL I. The sequencejjaf the peptide is 
SHRPPPPGHRV using the one-letter codesU. The four cen- 
tral prolin (?) residues form a rather rigid PP2 helix, the short 
peptide therefore does not show any folding or unfolding on 
the considered timescales. 

The atomic structure of the peptide was coarse-grained by 
first placing a van-der-Waals sphere around the Cq atoms of 
each of the residues. The side chains of the residues, which 
extended beyond this first sphere, were enclosed in a second 
(and third) sphere, for which the positions and radii were cho- 
sen to reproduce the van-der-Waals surface of the peptide as 
close as possible. In each residue the effective charges from 
the PDB structure with \q\ < 0.5 e were kept. Then, the Cq 
spheres were connected by springs. Thejdiffusion coefficients 
of the residues were taken from Ma etalU. Running BD sim- 
ulations of this first setup we determined the resulting mutual 
center-of-mass (CM) distances for each pair of residues and 
compared them to the respective distances derived from a 20 
ns MD simulation of the peptide in a water box. This refejis 
ence simulation was performed with the Gromacs packageU 
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FIG. 3: CM diffusion coefficient Dcm = ^^^J of the peptide vs. 
the length of the observation time interval AT. The arrow labeled 
with Trei indicates the typical velocity relaxation time of a single 
residue of ~ 40 fs. With LD the long-time value of Dcm is obtained 
only for AT > 100 Trei- With HI, Dcm is increased by a factor of 
2.7. 

using Ihe OPLS all-atom force fielcS and the TIP4P water 
modelU. The MD trajectory consequently includes the fi- 
nite damping of and the hydrodynamic coupling between the 
residues. As this first set of springs was not enough to re- 
produce all distance distributions correctly, additional springs 
between the CMs of the residues were added and their val- 
ues optimized manually until a sufficient agreement between 
the MD results and a standard BD simulation was achieved. 
With this parametrized bead-spring-model of the peptide we 
ran BD and LD simulations, each with and without hydro- 
dynamics. In all simulations, a single copy of the peptide was 
simulated in an unbounded box. The comparison of these four 
different simulation schemes shows the effects of each of the 
simplifications. 

The calculated relaxation times of the peptide residues 
range from Trei = 32 fs for the small glycin (G) to 50 fs 
for the larger arginine (R). For a reliable propagation we used 
a conservatively estimated small time step At = 10 fs in all 
simulations. 



A. Center-of-Mass Diffusion 

For the diffusion of the CM of the peptide two effects were 
observed. As for polymers the CM diffusion coefficient Dcm 
was increased when HI was included. This can be well seen in 
figure ^ which compares the x-component of the CM motion 
of the peptide obtained with the same sequence of random 
numbers for the four different propagation schemes. All four 
curves show the same trends, but the amplitudes from BDh-HI 
and LDh-HI are nearly twice as large as without HI. In this 
plot one can also see that the trajectories with LD are much 
smoother than with BD due to the finite relaxation time. For 
both BD with and without HI the power spectrum can be well 
fitted with a 1// behavior, where the amplitudes with HI are 
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At [ps] 

FIG. 4: Orientadonal relaxation of the end-to-end vector R^e as seen 
from the autocorrelation Cee from BD-l-HI (red points) and LD-l-HI 
simulations (black squares). The solid line indicates the fit to a 
stretched exponential as explained in the text. 

about 60% higher. With LD, however, the amplitudes are the 
same in the low-frequency range as in the corresponding BD 
simulations, but fluctuations that occur on timescale shorter 
than a picosecond are suppressed by about 30% in this one- 
dimensional representation (data not shown). 

Similar conclusions can be drawn from how the CM dif- 
fusion coefficient Dcm{^T) changes with the length of 
the observation time step AT. Figure ^ gives Dcm = 
(Ar^ (AT)) / GAT for the four different propagation schemes. 
Here, for the same trajectory, which was simulated at a 
timestep At 10 fs, the observation window AT was var- 
ied at which the CM displacements Ar(AT) were obtained. 
As expected, the BD simulations have no implicit time scale 
and, consequently, give the same Dcm for any AT. The 
creeping flow HI used here accelerates the diffusive motion 
without introducing a timescale of its own and Dcm is about 
2.7 times larger with hydrodynamics than without. For long 
observation intervals AT > 10 ps, the LD propagation re- 
produces the respective BD values for Dcm, whereas in the 
ballistic short time regime Dcm increases linearly with AT. 
The relaxation times of the individual residues are in the range 
Trei = 32 ... 50 fs, which is indicated in figure |[ However, it 
takes about two orders of magnitude longer, until Dcm from 
the LD propagation reaches the timescale-free BD value. 



B. Orientational Relaxation 

As shown above, the inclusion of HI accelerates the collec- 
tive CM diffusion. This is accompanied by a slowdown of the 
internal dynamics. For bead-spring polymers this can, e.g., be 
seen in a prolonged correlation of the end-to-end vector Ree 
pointing from the first to the last monomer The peptide used 
here is too short and not flexible enough to show pronounced 
differences in the relaxation of Ree whether HI are included 
or not. Nevertheless, the effects of the finite damping can be 
seen in the autocorrelation Cee(At) = {Ree{t)Ree{t + At)), 
which is plotted in figure Q as 1 — Cee (At). For both BD and 
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FIG. 5: Comparison of the power spectra for the dynamics of the 
CM distance between the first and the last residue with the different 
propagation schemes. In all four panels the black data (points) are 
from the MD simulation and the red (squares) from the respective 
coarse grained method. 

BDh-HI, Cee fits well to a stretched exponential 

Cee(At) -exp[-(At/T)"] (15) 

with r = 4.1 ns and a = 0.95. With the LD propagation, 
the dynamics on timescales shorter than a few picoseconds is 
slowed down both with and without HI. Correspondingly, in 
figure |4[ 1 — Cee from the Langevin simulation is smaller than 
with the overdamped Brownian propagation for short time in- 
tervals, while for longer intervals the same exponential decay 
is observed. 



C. Internal Dynamics 

The bonds of the coarse-grained representation of the pep- 
tide had been parametrized against the stationary distance dis- 
tributions without considering the actual dynamics, which are 
also influenced by the residue masses and the damping. Con- 
sequently, the dynamics will have different temporal signa- 
tures with the different propagation schemes. Figure || shows 
the power spectrum of the distance between the first and the 
last bond for the four integration schemes in comparison to 
the power spectrum extracted from the MD trajectory. These 
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already more stable than BD because of the damping at high 
frequencies. 

The respective runtimes for 100000 timesteps on a 2.8 GHz 
Pentium 4 CPU were 18.9 s for BD, 24.2 s for BD+Hl, 19.5 
s for LD, and 24.5 s for LD+HI. As explained above, keep- 
ing track of the velocity in the LD propagation only adds a 
minimal cost which already pays off when only the numeri- 
cal stability is taken into account and considerations about the 
regime of damping are ignored. In this small peptide with its 
35 effective charges and 38 bonds, the hydrodynamic interac- 
tions between the eleven residues only add a manageable 25% 
to the total runtime, whereby the direct interactions and the 
hydrodynamic coupling have the same 0{N'^) runtime scal- 
ing. 



FIG. 6: Variance of the distance distribution of the bond between 
residues 5 and 6 vs. the simulation timestep At for the four propa- 
gation schemes. 

power spectra are typical for the distances between residues 
located on opposite ends of the peptide. 

As seen in figure ||, the first four peaks between 0.002 ps^^ 
and 0.01 ps~^ are present in the LD result and can be recog- 
nized in the LDh-HI spectrum. BD and BDh-HI are damped 
too fast. Consequently, no real peak structure can be identi- 
fied above the strongly fluctuating background. For the other 
distances between residues on opposite ends of the peptide, a 
similar trend is found. 

These spectra now could be further improved by optimiz- 
ing the effective masses of the residues. However, already 
these results with the raw textbook-masses demonstrate that 
the short-time dynamics on the picosecond timescale are bet- 
ter reproduced with the finitely damped Langevin algorithm. 

D. Numerical Stability and Runtime 

The first test had been to compare the distance distribu- 
tions between the CMs of the residues for all four simulation 
schemes. As expected, for these static measures no differ- 
ences beyond the statistical uncertainties were found. Get- 
ting the same positions and widths with LD as with BD shows 
that the velocity relaxation occurs faster than one oscillation 
period in the interparticle potentials even for the rather hard 
springs that had to be used along the peptide backbone be- 
tween the Cq atoms. Otherwise one would expect a widening 
of the distance distributions and, in the extreme case, a bi- 
modal density as in a classical weakly damped oscillator. 

The stability of each integration scheme was then tested via 
the variances of the distance distributions for the hard springs, 
i.e., those distributions which have a small variance. Fig- 
ure ^ shows the representative behavior for the distance be- 
tween the fifth and the sixth residue, which are both part of 
the stiff PP2 helix in the central part of the peptide. With BD 
alone, timesteps up to 0.05 ps would have been okay, whereas 
the viscous damping due to HI allows for three times longer 
timesteps. The same trend is observed for LD, which itself is 



IV. SUMMARY AND CONCLUSIONS 

In this publication we investigated the effects of a finite 
damping and of hydrodynamic interactions in Brownian dy- 
namics simulations of a small biological peptide. 

Integrating Newton's equations of motion with a contin- 
uous viscous solvent over one timestep, we arrived at a 
Langevin type integration scheme which allows for the inclu- 
sion of our recently introduced efficient approximation for the 
hydrodynamic coupUngU. With a fast damping this Langevin 
scheme reduces to, the conventional BD algorithm of Ermak 
and McCammonliJ. Due to the analytic integration for one 
timestep the Langevin propagator is about as fast numerically 
as the simple BD integration, but does not have the concep- 
tual limitation that the timestep At has to be much longer 
than the velocity relaxation time Trei- This issue arises es- 
pecially in simulations of proteins of different sizes, because 
there At may be too close to Trei for the larger proteins, while 
for the smaller, faster proteins At is already too long for a 
stable propagation. With the LD propagator, the length of At 
is only limited by numerical accuracy and stabihty considera- 
tions. Actually, our test showed that only when At is at least 
two orders of magnitude longer than Trei, the finite damping 
may safely be neglected. For practical applications such a 
long integration timestep is nearly never feasible. 

Comparing the dynamics of a small peptide from our coarse 
grained BD and LD simulations to an atomistic molecular dy- 
namics trajectory, we found that the distance distributions be- 
tween the residues of the peptide can be reproduced with a 
bead-spring model with either propagation technique and the 
same force constants. However, when the dynamics of the 
relative motions are compared, LD with its more realistic fi- 
nite damping gives an overall better agreement of the internal 
dynamics. 

For the short and rather stiff peptide investigated here, the 
influence of hydrodynamic interactions on the internal mo- 
tion was small. The center-of-mass motion, however, became 
faster by nearly a factor of three. Compared to a plain BD sim- 
ulation, the slightly delayed orientational motion with the LD 
scheme together with the faster translation due to the HI may 
change the microscopic dynamics during, for example, the as- 
sociation of two proteins or the binding of a small peptide 
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to a larger protein. At close encounter the smoothed motion 
with LD may then further modify the obtained association and 
dissociation rates towards a more stable bound state. From a 
conceptual point of view we therefore recommend to include 
hydrodynamics and to use a Langevin propagation scheme for 
coarse grained protein simulations. 

From a numerical point of view, we find no difference in 
runtime for BD vs. LD and a moderate increase in runtime of 
about 30% for adding our approximate hydrodynamics with 
its O(iV^) runtime scaling. The coarse grained model of the 
peptide had about six times as many effective charges on the 
individual residues than there were effective spheres for the 
hydrodynamics. For peptide or protein models with a similar 
ratio of the numbers of the simple Coulomb compared to the 
more expensive hydrodynamic interactions, the runtime cost 
incurred by HI will be similar 

The finitely damped LD propagation was also more stable 
numerically, i.e., one could use longer integration timesteps 
than with the overdamped BD algorithm. The numerical sta- 
bility of both the BD and the LD schemes further increased 
with the viscous damping due to the HI. Actually, the longer 
timesteps that could be used with LD and HI more than com- 
pensated for the increased runtime. 

Here we used a simple peptide to demonstrate that coarse 
grained simulation methods need not stop at the level of Ein- 
stein's formulation of Brownian diffusion but that they can 
also be used very efficiently for much smaller systems, which 
are usually considered the realm of atomistic simulations or 
partially coarsened methods, where a few atoms each are 



summarized into one pseudo atomUU. The main difference, 
which makes fine-grained BD and LD simulations so much 
faster than the coarse grained versions of all-atom descrip- 
tions is that the very many solvent molecules do not have to 
be treated explicitly. 

Summarizing we could show that the LDh-HI simulation 
scheme, which includes more of the mechanically induced 
solvent effects than the conventional BD algorithm, can be 
formulated and implemented efficiently. Consequently, we 
are looking forward to applying the methods presented here to 
larger and more realistic systems in the near future. Potential 
applications will be the inve^gation of association funnels of 
peptide-jprotein encounterl— U, the folding of coarse grained 
proteinsU, and dense many-particle scenarios like the cytosol 
inside a biological cell. 

The LD and HI algorithms presented here have been imple- 
mented in a general purpose multi-particle BD/LD simulation 
library, which will be available from the corresponding author 
free of charge for academic use. 
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